scripts/Empirical example NatMath/NatMathOptimal.R

In joakimwallmark/PolyOptimalIRT:

# ----------- preliminaries -----------
library(TestGardener)

# ----------- read in data -----------
titlestr  <- "National Math Test"
U         <- scan("data/NatMath.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

#  data in score mode, convert to index mode
U = U + 1

#  treat this test as a rating scale
key     <- NULL

# Define the option score values for each item ----------------------------

scoreList <- vector("list",32) # option scores
scoreList[[ 1]] = c(0, 1)
scoreList[[ 2]] = c(0, 1)
scoreList[[ 3]] = c(0, 1)
scoreList[[ 4]] = c(0, 1)
scoreList[[ 5]] = c(0, 1)
scoreList[[ 6]] = c(0, 1, 2)
scoreList[[ 7]] = c(0, 1)
scoreList[[ 8]] = c(0, 1)
scoreList[[ 9]] = c(0, 1)
scoreList[[10]] = c(0, 1, 2)
scoreList[[11]] = c(0, 1)
scoreList[[12]] = c(0, 1)
scoreList[[13]] = c(0, 1)
scoreList[[14]] = c(0, 1, 2)
scoreList[[15]] = c(0, 1, 2)
scoreList[[16]] = c(0, 1)
scoreList[[17]] = c(0, 1, 2)
scoreList[[18]] = c(0, 1, 2)
scoreList[[19]] = c(0, 1, 2)
scoreList[[20]] = c(0, 1, 2)
scoreList[[21]] = c(0, 1, 2, 3)
scoreList[[22]] = c(0, 1, 2)
scoreList[[23]] = c(0, 1, 2, 3)
scoreList[[24]] = c(0, 1, 2)
scoreList[[25]] = c(0, 1, 2)
scoreList[[26]] = c(0, 1, 2)
scoreList[[27]] = c(0, 1, 2)
scoreList[[28]] = c(0, 1, 2)
scoreList[[29]] = c(0, 1, 2, 3)
scoreList[[30]] = c(0, 1, 2)
scoreList[[31]] = c(0, 1, 2)
scoreList[[32]] = c(0, 1, 2, 3, 4)

itemLab <- vector("list", 32)
itemLab[[ 1]] <-  'Question  1 '
itemLab[[ 2]] <-  'Question  2a'
itemLab[[ 3]] <-  'Question  2b'
itemLab[[ 4]] <-  'Question  3a'
itemLab[[ 5]] <-  'Question  3b'
itemLab[[ 6]] <-  'Question  4 '
itemLab[[ 7]] <-  'Question  5 '
itemLab[[ 8]] <-  'Question  6 '
itemLab[[ 9]] <-  'Question  7 '
itemLab[[10]] <-  'Question  8a'
itemLab[[11]] <-  'Question  8b'
itemLab[[12]] <-  'Question  9a'
itemLab[[13]] <-  'Question  9b'
itemLab[[14]] <-  'Question 10 '
itemLab[[15]] <-  'Question 11 '
itemLab[[16]] <-  'Question 12a'
itemLab[[17]] <-  'Question 12b'
itemLab[[18]] <-  'Question 13 '
itemLab[[19]] <-  'Question 14 '
itemLab[[20]] <-  'Question 15 '
itemLab[[21]] <-  'Question 16 '
itemLab[[22]] <-  'Question 17a'
itemLab[[23]] <-  'Question 17b'
itemLab[[24]] <-  'Question 18 '
itemLab[[25]] <-  'Question 19 '
itemLab[[26]] <-  'Question 20 '
itemLab[[27]] <-  'Question 21a'
itemLab[[28]] <-  'Question 21b'
itemLab[[29]] <-  'Question 22 '
itemLab[[30]] <-  'Question 23 '
itemLab[[31]] <-  'Question 24 '
itemLab[[32]] <-  'Question 25 '

optList <- list(itemLab=itemLab, optLab=NULL, optScr=scoreList)
maxScore <- sum(sapply(scoreList, max))

NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore))
NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore), NumBasis = 5, nbin = 12)
NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore), nbin = 16)
NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore), NumBasis = 5)
NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore), NumBasis = 9, nbin = 15)
NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore), NumBasis = 9)
NatMath_dataList <- make.dataList(U, key, optList, scrrng = c(0,maxScore), NumBasis = 7)

#  Set initial values that are required in the later analysis -----------------
#  compute the initial option surprisal curves using the
#  percentage ranks as initial estimates of theta
theta     <- NatMath_dataList$percntrnk
thetaQnt  <- NatMath_dataList$thetaQnt


#  Proceed through the cycles---------------------------------------------------
ncycle=10
AnalyzeResult <- Analyze(theta, thetaQnt, NatMath_dataList, ncycle=ncycle, itdisp=FALSE)
#AnalyzeResult.bas5.bin12 <- AnalyzeResult
AnalyzeResult.default <- AnalyzeResult
AnalyzeResult.basdef.bin12 <- AnalyzeResult
AnalyzeResult.bas9.bin12 <- AnalyzeResult
AnalyzeResult.bas5.bindef <- AnalyzeResult
AnalyzeResult.bas9.bindef <- AnalyzeResult
AnalyzeResult.bas7.bindef <- AnalyzeResult
#
# head(AnalyzeResult$parList[[1]]$WfdList[[20]]$Wbin)
# head(AnalyzeResult$parList[[1]]$WfdList[[20]]$Pbin)
# head(AnalyzeResult$parList[[1]]$binctr)
# length(AnalyzeResult$parList)

parList  <- AnalyzeResult$parList
meanHvec <- AnalyzeResult$meanHvec

icycle <- 10

NatMath_parListi  <- parList[[icycle]]

WfdList    <- NatMath_parListi$WfdList
Qvec       <- NatMath_parListi$Qvec
binctr     <- NatMath_parListi$binctr
theta      <- NatMath_parListi$theta
arclength  <- NatMath_parListi$arclength
alfine     <- NatMath_parListi$alfine

#  ----------------------------------------------------------------------------
#                   Plot surprisal curves for each test question
#  ----------------------------------------------------------------------------

#  plot both the probability and surprisal curves along with data points
Wbinsmth.plot(binctr, Qvec, WfdList, NatMath_dataList, Wrng=c(0,7), saveplot = F)

#save(NatMath_dataList, AnalyzeResult, file = "data/NatMath_fittedmodel.RData")